#load orginal version of GSEA:
source("https://rgsea.googlecode.com/svn/trunk/src/GSEA.1.0.R")

#load useful funcs of MC package:
source("https://mc-r.googlecode.com/svn/trunk/src/mcLite.R")

#load dependencies:
library.mc("Biobase","bioc")
library.mc("GSEABase","bioc")
library.mc("ArrayTools","bioc")
library.mc("hgu95av2.db","bioc")
library.mc("GO.db","bioc")
library.mc("illuminaHumanv3.db","bioc")

runGSEA<-function( df,                                      # Gene expression data frame
                   classes,                                 # Target classes
                   outdir=getwd(),                          # Output directory (default: getwd())
                   doc.string            = "GSEA",          # Documentation string used as a prefix to name result files (default: "GSEA.analysis")
                   non.interactive.run   = F,               # Run in interactive (i.e. R GUI) or batch (R command line) mode (default: F)
                   reshuffling.type      = "sample.labels", # Type of permutation reshuffling: "sample.labels" or "gene.labels" (default: "sample.labels" 
                   nperm                 = 1000,            # Number of random permutations (default: 1000)
                   weighted.score.type   =  1,              # Enrichment correlation-based weighting: 0=no weight (KS), 1= weigthed, 2 = over-weigthed (default: 1)
                   nom.p.val.threshold   = -1,              # Significance threshold for nominal p-vals for gene sets (default: -1, no thres)
                   fwer.p.val.threshold  = -1,              # Significance threshold for FWER p-vals for gene sets (default: -1, no thres)
                   fdr.q.val.threshold   = 0.25,            # Significance threshold for FDR q-vals for gene sets (default: 0.25)
                   topgs                 = 20,              # Besides those passing test, number of top scoring gene sets used for detailed reports (default: 10)
                   adjust.FDR.q.val      = F,               # Adjust the FDR q-vals (default: F)
                   gs.size.threshold.min = 15,              # Minimum size (in genes) for database gene sets to be considered (default: 25)
                   gs.size.threshold.max = 500,             # Maximum size (in genes) for database gene sets to be considered (default: 500)
                   reverse.sign          = F,               # Reverse direction of gene list (pos. enrichment becomes negative, etc.) (default: F)
                   preproc.type          = 0,               # Preproc.normalization: 0=none, 1=col(z-score)., 2=col(rank) and row(z-score)., 3=col(rank). (def: 0)
                   random.seed           = 3338,            # Random number generator seed. (default: 123456)
                   perm.type             = 0,               # For experts only. Permutation type: 0 = unbalanced, 1 = balanced (default: 0)
                   fraction              = 1.0,             # For experts only. Subsampling fraction. Set to 1.0 (no resampling) (default: 1.0)
                   replace               = F,               # For experts only, Resampling mode (replacement or not replacement) (default: F)
                   save.intermediate.results = F,           # For experts only, save intermediate results (e.g. matrix of random perm. scores) (default: F)
                   OLD.GSEA              = F,               # Use original (old) version of GSEA (default: F)
                   use.fast.enrichment.routine = T          # Use faster routine to compute enrichment for random permutations (default: T)
                   ){
  
  print("[rGSEA] Creating GSEA files...")
  eSet <- ExpressionSet(assayData=df)
  eSet$Class<-classes
  createGSEAFiles(mydir=outdir,eSet, "Class")
  
  print("[rGSEA] Creating gene set collections...")


  if (grepl("at",rownames(df)[1])){#For Affy Probe IDs
    gsc <- GeneSetCollection(rownames(df),
                             idType=AnnotationIdentifier("hgu95av2"),
                             setType=GOCollection())
  } else if (grepl("ILMN",rownames(df)[1])) {#For Illumina Probe IDs
    gsc <- GeneSetCollection(rownames(df),
                           idType=AnnotationIdentifier("illuminaHumanv3"),
                           setType=GOCollection())
  } else {#For ENTREZ IDs
    gsc <- GeneSetCollection(rownames(df),
                           idType=AnnotationIdentifier("org.Hs.eg.db"),
                           setType=GOCollection())
  }


  print("[rGSEA] Creating gene set database...")
  fl <- paste(outdir,"/Result/Pathway_Analysis/GSEA/GS.db.gmt",sep="")
  toGmt(gsc, fl)
  

    
  GSEA(                        
    input.ds = paste(outdir,"/Result/Pathway_Analysis/GSEA/Class.probe.gct",sep=""), 
    input.cls = paste(outdir,"/Result/Pathway_Analysis/GSEA/Class.phenotype.cls",sep=""), 
    gs.db = paste(outdir,"/Result/Pathway_Analysis/GSEA/GS.db.gmt",sep=""), 
    output.directory = paste(outdir,"/Result/Pathway_Analysis/GSEA/",sep=""), #needs trailing /!
     #  Program parameters :-------------------------------------------------------------------------------------------------------------------------
     doc.string            = doc.string,   
     non.interactive.run   = non.interactive.run,               
     reshuffling.type      = reshuffling.type, 
     nperm                 = nperm,            
     weighted.score.type   =  weighted.score.type,              
     nom.p.val.threshold   = nom.p.val.threshold,              
     fwer.p.val.threshold  = fwer.p.val.threshold,              
     fdr.q.val.threshold   = fdr.q.val.threshold,            
     topgs                 = topgs,              
     adjust.FDR.q.val      = adjust.FDR.q.val,               
     gs.size.threshold.min = gs.size.threshold.min,              
     gs.size.threshold.max = gs.size.threshold.max,             
     reverse.sign          = reverse.sign,               
     preproc.type          = preproc.type,               
     random.seed           = random.seed,            
     perm.type             = perm.type,               
     fraction              = fraction,             
     replace               = replace,               
     save.intermediate.results = save.intermediate.results,           
     OLD.GSEA              = OLD.GSEA,               
     use.fast.enrichment.routine = use.fast.enrichment.routine 
  )

 
  
  GSEA.Analyze.Sets( 
    directory =  paste(outdir,"/Result/Pathway_Analysis/GSEA/",sep=""), #need trailing / 
    topgs = 20, 
    height = 16, 
    width = 16 
  ) 
  
 
  
}

  

